################################################################################################################
#                                         Reading data                                                          
################################################################################################################
TAB.IEA <- read.csv("D:/Arbeit/Forschung_aktuell/Environmental_Treaties_Paper/Data_Analysis/Replication-files/TAB.IEA.csv")
str(TAB.IEA, list.len = ncol(TAB.IEA)) 
                                       

####################################################################################################################
#              Deselect treaties where there was no singanture stage with none_signed==0                            
####################################################################################################################
TAB.IEA <- subset(TAB.IEA, TAB.IEA$none_signed==0)


################ experience_i_sign ##########################################
    TAB.IEA$summed.exp.2.std <- TAB.IEA$summed.exp.2
    TAB.IEA$summed.exp.2.std <- (TAB.IEA$summed.exp.2.std - mean(TAB.IEA$summed.exp.2.std, na.rm=TRUE)) / sd(TAB.IEA$summed.exp.2.std, na.rm=TRUE)
##########################################################################

################ raw number of countries' signatures ##########################################
tmp <- by(TAB.IEA, TAB.IEA$treaty, function (x) length(unique(x$cowcode[x$first_signer == 1])))
TAB.IEA$n_i_sign <- tmp[as.character(TAB.IEA$treaty)] 


####################################################################################################################
#              Deselecting protocols & amendments                                                                   
####################################################################################################################
TAB.IEA.nopro <- subset(TAB.IEA, TAB.IEA$prot.amd==0)
unique(TAB.IEA.nopro$treaty)

########################################################################################################################
#                                            Descriptive statistics                                                     
########################################################################################################################

   summary(TAB.IEA.nopro$docs_i_sign)
   summary(TAB.IEA.nopro$patent_i_sign)
   summary(TAB.IEA.nopro$summed.exp.2.std)
   
   summary(TAB.IEA.nopro$first_signer)
   summary(TAB.IEA.nopro$wealth_i_sign)
   summary(TAB.IEA.nopro$power_i_sign)
   summary(TAB.IEA.nopro$r_coast_land)         
   summary(TAB.IEA.nopro$IO_membership)
   summary(TAB.IEA.nopro$threshold)                                                                 
   summary(TAB.IEA.nopro$lagpercregion)   
   summary(TAB.IEA.nopro$open)       
   summary(TAB.IEA.nopro$rgdpl)      
   summary(TAB.IEA.nopro$rgdplsq)  
   summary(TAB.IEA.nopro$lnso2pc)               
   summary(TAB.IEA.nopro$meanpc)   
   summary(TAB.IEA.nopro$polity2) 
   summary(TAB.IEA.nopro$GDPL)     


   length(TAB.IEA.nopro$docs_i_sign[is.na(TAB.IEA.nopro$docs_i_sign)==FALSE])
   length(TAB.IEA.nopro$patent_i_sign[is.na(TAB.IEA.nopro$patent_i_sign)==FALSE])  
   length(TAB.IEA.nopro$summed.exp.2.std[is.na(TAB.IEA.nopro$summed.exp.2.std)==FALSE])
   
   length(TAB.IEA.nopro$first_signer[is.na(TAB.IEA.nopro$first_signer)==FALSE])   
   length(TAB.IEA.nopro$wealth_i_sign[is.na(TAB.IEA.nopro$wealth_i_sign)==FALSE])
   length(TAB.IEA.nopro$power_i_sign[is.na(TAB.IEA.nopro$power_i_sign)==FALSE])
   length(TAB.IEA.nopro$r_coast_land[is.na(TAB.IEA.nopro$r_coast_land)==FALSE])
   length(TAB.IEA.nopro$IO_membership[is.na(TAB.IEA.nopro$IO_membership)==FALSE]) 
   length(TAB.IEA.nopro$threshold[is.na(TAB.IEA.nopro$threshold)==FALSE])   
   length(TAB.IEA.nopro$lagpercregion[is.na(TAB.IEA.nopro$lagpercregion)==FALSE])
   length(TAB.IEA.nopro$open[is.na(TAB.IEA.nopro$open)==FALSE])
   length(TAB.IEA.nopro$rgdpl[is.na(TAB.IEA.nopro$rgdpl)==FALSE])
   length(TAB.IEA.nopro$rgdplsq[is.na(TAB.IEA.nopro$rgdplsq)==FALSE])
   length(TAB.IEA.nopro$lnso2pc[is.na(TAB.IEA.nopro$lnso2pc)==FALSE])
   length(TAB.IEA.nopro$meanpc[is.na(TAB.IEA.nopro$meanpc)==FALSE])
   length(TAB.IEA.nopro$polity2[is.na(TAB.IEA.nopro$polity2)==FALSE])  
   length(TAB.IEA.nopro$GDPL[is.na(TAB.IEA.nopro$GDPL)==FALSE])


   sd(TAB.IEA.nopro$docs_i_sign, na.rm=TRUE)    
   sd(TAB.IEA.nopro$patent_i_sign, na.rm=TRUE)    
   sd(TAB.IEA.nopro$summed.exp.2.std, na.rm=TRUE)
   
   sd(TAB.IEA.nopro$first_signer, na.rm=TRUE)
   sd(TAB.IEA.nopro$wealth_i_sign, na.rm=TRUE)
   sd(TAB.IEA.nopro$power_i_sign, na.rm=TRUE)
   sd(TAB.IEA.nopro$r_coast_land, na.rm=TRUE)
   sd(TAB.IEA.nopro$IO_membership, na.rm=TRUE) 
   sd(TAB.IEA.nopro$threshold, na.rm=TRUE)
   sd(TAB.IEA.nopro$lagpercregion, na.rm=TRUE)
   sd(TAB.IEA.nopro$open, na.rm=TRUE)  
   sd(TAB.IEA.nopro$rgdpl, na.rm=TRUE)
   sd(TAB.IEA.nopro$rgdplsq, na.rm=TRUE)
   sd(TAB.IEA.nopro$lnso2pc, na.rm=TRUE)
   sd(TAB.IEA.nopro$meanpc, na.rm=TRUE)  
   sd(TAB.IEA.nopro$GDPL, na.rm=TRUE)
   sd(TAB.IEA.nopro$polity2, na.rm=TRUE)

##########################################################################################################################
#                                             Correlation matrix                                                           
##########################################################################################################################
TAB.temp <- subset(TAB.IEA.nopro, select=c("cowcode","stop","first_signer","docs_i_sign","patent_i_sign","summed.exp.2.std","marine_i_signer",
                                            "wealth_i_sign","power_i_sign"
                                           ))                                               
TAB.temp <- data.frame(unique(cbind(TAB.temp$cowcode,
                                    TAB.temp$stop,
                                    TAB.temp$first_signer,
                                    TAB.temp$docs_i_sign,
                                    TAB.temp$patent_i_sign,
                                    TAB.temp$summed.exp.2.std,
                                    TAB.temp$marine_i_signer,
                                    TAB.temp$wealth_i_sign,
                                    TAB.temp$power_i_sign
                                    )))
                                    
colnames(TAB.temp)[1:ncol(TAB.temp)] <- c("cowcode","stop","first_signer","docs_i_sign","patent_i_sign","summed.exp.2.std","marine_i_signer",
                                          "wealth_i_sign","power_i_sign"
                                          )                                                

for(i in 2:ncol(TAB.temp)){
    TAB.temp[,i] <- as.numeric(as.character(TAB.temp[,i]))
}
pairs(TAB.temp[,3:ncol(TAB.temp)])

temp <- subset(TAB.temp, 
            is.na(TAB.temp$first_signer)==FALSE & 
            is.na(TAB.temp$docs_i_sign)==FALSE &
            is.na(TAB.temp$patent_i_sign)==FALSE & 
            is.na(TAB.temp$summed.exp.2.std)==FALSE &
            is.na(TAB.temp$marine_i_sign)==FALSE &
            is.na(TAB.temp$wealth_i_sign)==FALSE & 
            is.na(TAB.temp$power_i_sign)==FALSE,
select=c("first_signer","docs_i_sign","patent_i_sign","summed.exp.2.std","marine_i_signer","wealth_i_sign","power_i_sign"
))
cor(temp)   
   
##################################################################################################################
#                                        Survival fits                                                            
##################################################################################################################
library(survival)

    Fit.1 <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                                first_signer
                                + patent_i_sign
                                        + IO_membership
                                        + threshold
                                        + lagpercregion
                                            + open
                                            + rgdpl
                                            + rgdplsq
                                            + lnso2pc
                                            + meanpc
                                            + GDPL 
                                            
                                + cluster(cowcode),data=TAB.IEA.nopro)
                summary(Fit.1)
                # no. states
                temp <-   subset(TAB.IEA.nopro,
                                is.na(first_signer)==FALSE &
                                is.na(patent_i_sign)==FALSE &
                                    is.na(IO_membership)==FALSE &   
                                    is.na(threshold)==FALSE &
                                    is.na(lagpercregion)==FALSE &
                                        is.na(open)==FALSE &
                                        is.na(rgdpl)==FALSE &
                                        is.na(rgdplsq)==FALSE &
                                        is.na(lnso2pc)==FALSE &
                                        is.na(meanpc)==FALSE &
                                        is.na(GDPL)==FALSE 
                            )               
                length(unique(temp$cowcode))
                # period of analysis
                min(temp$stop, na.rm=TRUE)    
                max(temp$stop, na.rm=TRUE)

# adding number of countries' signatures
(Fit.1.ns <- update(Fit.1, . ~ . + n_i_sign))
# number of countries' signatures, interacted with main independent variable
(Fit.1.int <- update(Fit.1, . ~ . + n_i_sign*patent_i_sign))
          
            
    Fit.2 <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                                first_signer 
                                + summed.exp.2.std
                                        + IO_membership
                                        + threshold
                                        + lagpercregion
                                            + open
                                            + rgdpl
                                            + rgdplsq
                                            + lnso2pc
                                            + meanpc
                                            + GDPL 
                                + cluster(cowcode),data=TAB.IEA.nopro)
                summary(Fit.2)
                    # no. states
                temp <-   subset(TAB.IEA.nopro,
                                is.na(summed.exp.2.std)==FALSE &
                                is.na(first_signer)==FALSE &
                                    is.na(IO_membership)==FALSE &   
                                    is.na(threshold)==FALSE &
                                    is.na(lagpercregion)==FALSE &
                                        is.na(open)==FALSE &
                                        is.na(rgdpl)==FALSE &
                                        is.na(rgdplsq)==FALSE &
                                        is.na(lnso2pc)==FALSE &
                                        is.na(meanpc)==FALSE &
                                        is.na(GDPL)==FALSE   
                            )               
                length(unique(temp$cowcode))
                # period of analysis
                min(temp$stop, na.rm=TRUE)    
                max(temp$stop, na.rm=TRUE)  

# number of countries' signatures
(Fit.2.ns <- update(Fit.2, . ~ . + n_i_sign))
# number of countries' signatures, interacted with main independent variable
(Fit.2.int <- update(Fit.2, . ~ . + n_i_sign*summed.exp.2.std))  
    
    Fit.4 <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                                first_signer
                                + patent_i_sign 
                                + wealth_i_sign
                                + power_i_sign 
                                        + IO_membership
                                        + threshold
                                        + lagpercregion
                                            + open
                                            + rgdpl
                                            + rgdplsq
                                            + lnso2pc
                                            + meanpc
                                            + GDPL 
                                + cluster(cowcode),data=TAB.IEA.nopro)
                summary(Fit.4)
                # no. states
                temp <-   subset(TAB.IEA.nopro,
                                is.na(first_signer)==FALSE &
                                is.na(patent_i_sign)==FALSE &
                                is.na(wealth_i_sign)==FALSE &
                                is.na(power_i_sign)==FALSE &
                                    is.na(IO_membership)==FALSE &   
                                    is.na(threshold)==FALSE &
                                    is.na(lagpercregion)==FALSE &
                                        is.na(open)==FALSE &
                                        is.na(rgdpl)==FALSE &
                                        is.na(rgdplsq)==FALSE &
                                        is.na(lnso2pc)==FALSE &
                                        is.na(meanpc)==FALSE &
                                        is.na(GDPL)==FALSE 
                            )               
                length(unique(temp$cowcode))
                # period of analysis
                min(temp$stop, na.rm=TRUE)    
                max(temp$stop, na.rm=TRUE)     

# adding number of countries' signatures
(Fit.4.ns <- update(Fit.4, . ~ . + n_i_sign))
# number of countries' signatures, interacted with main independent variable
(Fit.4.int <- update(Fit.4, . ~ . + n_i_sign*patent_i_sign))
                        
  
            ###  Marine model
            TAB.mar.fish.nopro <- subset(TAB.IEA.nopro, TAB.IEA.nopro$mar_fish==1)
            
            length(TAB.mar.fish.nopro$marine_i_signer[is.na(TAB.mar.fish.nopro$marine_i_signer)==FALSE])
            length(TAB.mar.fish.nopro$r_coast_land[is.na(TAB.mar.fish.nopro$r_coast_land)==FALSE])
            summary(TAB.mar.fish.nopro$marine_i_signer, digits = 3)
            summary(TAB.mar.fish.nopro$r_coast_land, digits = 3)
            sd(TAB.mar.fish.nopro$marine_i_signer, na.rm=TRUE) 
            sd(TAB.mar.fish.nopro$r_coast_land, na.rm=TRUE)
            
    Fit.3 <- coxph( Surv(TAB.mar.fish.nopro$start,TAB.mar.fish.nopro$stop,TAB.mar.fish.nopro$status)~
                                first_signer
                                + marine_i_signer
                                + r_coast_land 
                                        + IO_membership
                                        + threshold
                                        + lagpercregion
                                            + open
                                            + rgdpl
                                            + rgdplsq
                                            + lnso2pc
                                            + meanpc
                                            + GDPL                        
                                + cluster(cowcode),data=TAB.mar.fish.nopro)
                summary(Fit.3)
                # no. states
                temp <-   subset(TAB.mar.fish.nopro,
                                is.na(first_signer)==FALSE &
                                is.na(marine_i_signer)==FALSE &
                                is.na(r_coast_land)==FALSE &
                                    is.na(IO_membership)==FALSE &   
                                    is.na(threshold)==FALSE &
                                    is.na(lagpercregion)==FALSE &
                                        is.na(open)==FALSE &
                                        is.na(rgdpl)==FALSE &
                                        is.na(rgdplsq)==FALSE &
                                        is.na(lnso2pc)==FALSE &
                                        is.na(meanpc)==FALSE &
                                        is.na(GDPL)==FALSE   
                                        )               
                length(unique(temp$cowcode))
                # period of analysis
                min(temp$stop, na.rm=TRUE)    
                max(temp$stop, na.rm=TRUE) 
 
# adding number of countries' signatures
(Fit.3.ns <- update(Fit.3, . ~ . + n_i_sign))
# number of countries' signatures, interacted with main independent variable
(Fit.3.int <- update(Fit.3, . ~ . + n_i_sign*marine_i_signer))
 
                
    
    ### A1 - Patent model without controls
    Fit.1_0 <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                            first_signer
                            + patent_i_sign  
                            + cluster(cowcode),data=TAB.IEA.nopro)
            summary(Fit.1_0)
            # no. states
            temp <-   subset(TAB.IEA.nopro,
                            is.na(first_signer)==FALSE &
                            is.na(patent_i_sign)==FALSE 
                        )               
            length(unique(temp$cowcode))
            # period of analysis
            min(temp$stop, na.rm=TRUE)    
            max(temp$stop, na.rm=TRUE) 
      
        
    ### A2 Experience model without controls
    Fit.2_0 <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                            first_signer 
                            + summed.exp.2.std
                            + cluster(cowcode),data=TAB.IEA.nopro)
            summary(Fit.2_0)
                # no. states
            temp <-   subset(TAB.IEA.nopro,
                            is.na(summed.exp.2.std)==FALSE &
                            is.na(first_signer)==FALSE   
                        )               
            length(unique(temp$cowcode))
            # period of analysis
            min(temp$stop, na.rm=TRUE)    
            max(temp$stop, na.rm=TRUE)    
     
            
    ### A3 Marine model without controls
    Fit.3_0 <- coxph( Surv(TAB.mar.fish.nopro$start,TAB.mar.fish.nopro$stop,TAB.mar.fish.nopro$status)~
                            first_signer
                            + marine_i_signer
                            + r_coast_land 
                            + cluster(cowcode),data=TAB.mar.fish.nopro)
            summary(Fit.3_0)
            # no. states
            temp <-   subset(TAB.mar.fish.nopro,
                            is.na(first_signer)==FALSE &
                            is.na(marine_i_signer)==FALSE &
                            is.na(r_coast_land)==FALSE  
                                    )               
            length(unique(temp$cowcode))
            # period of analysis
            min(temp$stop, na.rm=TRUE)    
            max(temp$stop, na.rm=TRUE) 
       
            
    ### A4 Alternative model without controls
    Fit.4_0 <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                            first_signer
                            + patent_i_sign 
                            + wealth_i_sign
                            + power_i_sign 
                            + cluster(cowcode),data=TAB.IEA.nopro)
            summary(Fit.4_0)
                # no. states
            temp <-   subset(TAB.IEA.nopro,
                            is.na(first_signer)==FALSE &
                            is.na(patent_i_sign)==FALSE &
                            is.na(wealth_i_sign)==FALSE &
                            is.na(power_i_sign)==FALSE 
                        )               
            length(unique(temp$cowcode))
            # period of analysis
            min(temp$stop, na.rm=TRUE)    
            max(temp$stop, na.rm=TRUE) 
        
    ### A5 stratified on issues
    ############################################################################################
                    #1   General/Governance
                    #2   Atmosphere/Air/Outer Space
                    #3   Hazardous substances
                    #4   Marine environment/pollution
                    #5   Sustainabilty/Nature conservation and terrestrial living resources
                    #6   Energy
                    #7   Nuclear safety
                    #8   Marine living resources
                    #9   Freshwater sources
                    #10  Conflict/disasters
    ############################################################################################
    Fit.1.strat <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                                first_signer
                                + patent_i_sign 
                                        + IO_membership
                                        + threshold
                                        + lagpercregion
                                            + open
                                            + rgdpl
                                            + rgdplsq
                                            + lnso2pc
                                            + meanpc
                                            + GDPL 
                                + strata(issues)
                                + cluster(cowcode),data=TAB.IEA.nopro)
                summary(Fit.1.strat)
                # no. states
                temp <-   subset(TAB.IEA.nopro,
                                is.na(first_signer)==FALSE &
                                is.na(patent_i_sign)==FALSE &
                                    is.na(IO_membership)==FALSE &   
                                    is.na(threshold)==FALSE &
                                    is.na(lagpercregion)==FALSE &
                                        is.na(open)==FALSE &
                                        is.na(rgdpl)==FALSE &
                                        is.na(rgdplsq)==FALSE &
                                        is.na(lnso2pc)==FALSE &
                                        is.na(meanpc)==FALSE &
                                        is.na(GDPL)==FALSE &
                                        is.na(issues)==FALSE
                            )               
                length(unique(temp$cowcode))
                # period of analysis
                min(temp$stop, na.rm=TRUE)    
                max(temp$stop, na.rm=TRUE)  
     
     
    ### A6 clustered on 'treaty' instead of 'country'              
    Fit.1.treat.cl <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                            first_signer
                            + patent_i_sign 
                                    + IO_membership
                                    + threshold
                                    + lagpercregion
                                        + open
                                        + rgdpl
                                        + rgdplsq
                                        + lnso2pc
                                        + meanpc
                                        + GDPL 
                            + cluster(treaty),data=TAB.IEA.nopro)
            summary(Fit.1.treat.cl)
            # no. states
            temp <-   subset(TAB.IEA.nopro,
                            is.na(first_signer)==FALSE &
                            is.na(patent_i_sign)==FALSE &
                                is.na(IO_membership)==FALSE &   
                                is.na(threshold)==FALSE &
                                is.na(lagpercregion)==FALSE &
                                    is.na(open)==FALSE &
                                    is.na(rgdpl)==FALSE &
                                    is.na(rgdplsq)==FALSE &
                                    is.na(lnso2pc)==FALSE &
                                    is.na(meanpc)==FALSE &
                                    is.na(GDPL)==FALSE 
                        )               
            length(unique(temp$treaty))
            # period of analysis
            min(temp$stop, na.rm=TRUE)    
            max(temp$stop, na.rm=TRUE)  
     
     
    ### A7 including polity2 instead of meanpc 
    Fit.1.pol <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                       first_signer
                       + patent_i_sign 
                            + IO_membership
                            + threshold
                            + lagpercregion
                                + open
                                + rgdpl
                                + rgdplsq
                                + lnso2pc
                                + polity2
                                + GDPL 
                       + cluster(cowcode),data=TAB.IEA.nopro)
                summary(Fit.1.pol)
                # no. states
                temp <-   subset(TAB.IEA.nopro,
                                is.na(first_signer)==FALSE &
                                is.na(patent_i_sign)==FALSE &
                                    is.na(IO_membership)==FALSE &   
                                    is.na(threshold)==FALSE &
                                    is.na(lagpercregion)==FALSE &
                                        is.na(open)==FALSE &
                                        is.na(rgdpl)==FALSE &
                                        is.na(rgdplsq)==FALSE &
                                        is.na(lnso2pc)==FALSE &
                                        is.na(polity2)==FALSE &
                                        is.na(GDPL)==FALSE 
                            )               
                length(unique(temp$cowcode))
                # period of analysis
                min(temp$stop, na.rm=TRUE)    
                max(temp$stop, na.rm=TRUE)   
        
    ### A8 without first_signer
    Fit.1.wo.fs <- coxph( Surv(TAB.IEA.nopro$start,TAB.IEA.nopro$stop,TAB.IEA.nopro$status)~
                            patent_i_sign 
                                    + IO_membership
                                    + threshold
                                    + lagpercregion
                                        + open
                                        + rgdpl
                                        + rgdplsq
                                        + lnso2pc
                                        + meanpc
                                        + GDPL 
                            + cluster(cowcode),data=TAB.IEA.nopro)
            summary(Fit.1.wo.fs)
            # no. states
            temp <-   subset(TAB.IEA.nopro,
                            is.na(patent_i_sign)==FALSE &
                                is.na(IO_membership)==FALSE &   
                                is.na(threshold)==FALSE &
                                is.na(lagpercregion)==FALSE &
                                    is.na(open)==FALSE &
                                    is.na(rgdpl)==FALSE &
                                    is.na(rgdplsq)==FALSE &
                                    is.na(lnso2pc)==FALSE &
                                    is.na(meanpc)==FALSE &
                                    is.na(GDPL)==FALSE 
                        )               
            length(unique(temp$cowcode))
            # period of analysis
            min(temp$stop, na.rm=TRUE)    
            max(temp$stop, na.rm=TRUE) 
        
###################################################################################################################################
#                                    Simulation                                                                                    
###################################################################################################################################

N.Boot <- 199
Boot.i<-1
for(Boot.i in 0:N.Boot){
    
    start.time <- Sys.time()
    
    if(Boot.i>0){

        treaty.vec <- sort(unique(TAB.IEA.nopro$treaty))
        treaty.sample <- sample(1:length(treaty.vec), size=length(treaty.vec), replace = TRUE)
        treaty.sample.vec <- treaty.vec[treaty.sample]
        
        TAB.sim <- data.frame(treaty=treaty.sample.vec)
        TAB.sim <- aggregate(rep(1, nrow(TAB.sim)), by=list(TAB.sim$treaty), FUN=sum)
        colnames(TAB.sim)[1:2] <- c("treaty","n")
        
        n.vec <- sort(unique(TAB.sim$n))
     
        for(n.i in n.vec){
        
            temp <- subset(TAB.sim, TAB.sim$n==n.i, select=-n)
            temp$bootselect <- 1
            
            for(n.j in 1:n.i){
            
                TAB.IEA.boot.temp <- merge(TAB.IEA.nopro, temp, all.x=TRUE, all.y=FALSE, by="treaty")
                TAB.IEA.boot.temp <- subset(TAB.IEA.boot.temp, TAB.IEA.boot.temp$bootselect==1, select=-bootselect)
                
                if(n.i==n.vec[1] & n.j==n.i[1]){
                    TAB.IEA.boot <- TAB.IEA.boot.temp
                }else{
                    TAB.IEA.boot <- rbind(TAB.IEA.boot, TAB.IEA.boot.temp)
                }
            
            }
            
        }
        
     }else{
     
        TAB.IEA.boot <- TAB.IEA.nopro
     
     } 
        
    if((class(try(
 
        Fit <- coxph( Surv(start,stop,status)~
                          open
                        + meanpc
                        + IO_membership
                        + rgdpl
                        + rgdplsq
                        + lnso2pc
                        + threshold
                        + lagpercregion
                        + GDPL
                        + first_signer
                        + patent_i_sign    
                        + cluster(cowcode),data=TAB.IEA.boot)
                                        
    ))=="try-error")==FALSE){  
    
        temp.coef <- data.frame(summary(Fit)$coefficients)    
        temp.coef$boot <- Boot.i
        temp.coef$varname <- rownames(temp.coef)
        rownames(temp.coef) <- NULL
        
        if(Boot.i==0){
            TAB.coef <- temp.coef
        }else{
            TAB.coef <- rbind(TAB.coef, temp.coef)
        }
        
        end.time <- Sys.time()
        duration.time <- end.time - start.time
        
        cat("Simulation ", Boot.i, " of ", N.Boot, " |  Duration ", duration.time, "\n")
        
    }
    
    write.csv(TAB.coef, file="D:/Arbeit/Forschung_aktuell/Environmental_Treaties_Paper/Data_Analysis/Replication-filesc/TAB_Simulation_Fit.1.csv", row.names=FALSE)

} 

TAB.Result <- read.csv(file="D:/Arbeit/Forschung_aktuell/Environmental_Treaties_Paper/Data_Analysis/Replication-files/TAB_Simulation_Fit.1.csv")

Var.Nam.vec <- c("first_signer","patent_i_sign") 

i <- 2

for(i in 1:length(Var.Nam.vec)){
    
    Var.Nam.i <- Var.Nam.vec[i]
    coef.est <- TAB.Result$coef[TAB.Result$varname==Var.Nam.i & TAB.Result$boot==0]    
    coef.values <- TAB.Result$coef[TAB.Result$varname==Var.Nam.i & TAB.Result$boot!=0 & is.na(TAB.Result$coef)==FALSE]
    coef.mean <- mean(coef.values)
    coef.max <- max(coef.values)    
    coef.min <- min(coef.values)            
    coef.median <- median(coef.values)
    coef.P.0.05 <- quantile(coef.values, probs = 0.05, na.rm = TRUE)
    coef.P.0.95 <- quantile(coef.values, probs = 0.95, na.rm = TRUE)
    coef.P.0.025 <- quantile(coef.values, probs = 0.025, na.rm = TRUE)
    coef.P.0.975 <- quantile(coef.values, probs = 0.975, na.rm = TRUE)
    coef.robust.se <- TAB.Result$robust.se[TAB.Result$varname==Var.Nam.i & TAB.Result$boot==0]
    coef.ci.95.up <- coef.est + qt(p=(1-0.05/2), df=Inf) * coef.robust.se
    coef.ci.95.lo <- coef.est - qt(p=(1-0.05/2), df=Inf) * coef.robust.se    
    
    coef.se <- TAB.Result$se.coef.[TAB.Result$varname==Var.Nam.i & TAB.Result$boot==0]
    coef.sd <- sd(coef.values)
    
    plot.range <- max(abs(c(coef.max,coef.min)))
    
    TAB.density <- density(coef.values)

    y.max <- max(TAB.density$y, na.rm=TRUE)
    
    par(mfrow=c(1,1), mar=c(3,0.5,3,0.5)
        , family="mono", cex=1.5, yaxt="n", xaxt="s")

    plot(NULL, NULL, xlim=c(-plot.range,plot.range), ylim=c(0,y.max),  main=Var.Nam.i, xlab="", ylab="")
    
    points(TAB.density$x, TAB.density$y, type="l", lwd=2, lty=1)
    
    abline(v=coef.mean, lty=1, lwd=1, col=gray(0.6))
    abline(v=coef.median, lty=2, lwd=1, col=gray(0.6))
    abline(v=coef.est,  lty=1, lwd=1, col=1)       
    
    abline(v=coef.P.0.025,  lty=2, lwd=1, col=1)
    abline(v=coef.P.0.975,  lty=2, lwd=1, col=1)
    
    legend("topleft",
           legend=c("estimate","mean","median","95% intv."), 
           lty=c(1,1,2,2),
           lwd=c(1,1,1,1), 
           col=c(1,gray(0.6),gray(0.6),1),
           bty = "n",            
           cex = 1                                  
    )    
        
    dev.print(device=pdf, file=paste("D:/Arbeit/Forschung_aktuell/Environmental_Treaties_Paper/Data_Analysis/Replication-files/",Var.Nam.i ,".pdf",sep=""))                    
   
}
